Cargo librerias
library(spatstat.core)
## Loading required package: spatstat.data
## Loading required package: spatstat.geom
## spatstat.geom 2.3-0
## Loading required package: nlme
## Loading required package: rpart
## spatstat.core 2.3-1
library(plotly)
## Loading required package: ggplot2
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
Principio de Aceptacion-Rechazo para un Proceso Binomial
# Graficacion
cuadri <- function(x, y){sqrt(x^2+y^2)/sqrt(8)} # Funcion de intensidad, toma valores entre 0 y 1
n<-10 # cant de datos
x <- y <- seq(0, 2, length= n)
xx<-rep(x,n)
yy<-rep(y,rep(n,n))
cbind(xx,yy) # grilla en el soporte
## xx yy
## [1,] 0.0000000 0.0000000
## [2,] 0.2222222 0.0000000
## [3,] 0.4444444 0.0000000
## [4,] 0.6666667 0.0000000
## [5,] 0.8888889 0.0000000
## [6,] 1.1111111 0.0000000
## [7,] 1.3333333 0.0000000
## [8,] 1.5555556 0.0000000
## [9,] 1.7777778 0.0000000
## [10,] 2.0000000 0.0000000
## [11,] 0.0000000 0.2222222
## [12,] 0.2222222 0.2222222
## [13,] 0.4444444 0.2222222
## [14,] 0.6666667 0.2222222
## [15,] 0.8888889 0.2222222
## [16,] 1.1111111 0.2222222
## [17,] 1.3333333 0.2222222
## [18,] 1.5555556 0.2222222
## [19,] 1.7777778 0.2222222
## [20,] 2.0000000 0.2222222
## [21,] 0.0000000 0.4444444
## [22,] 0.2222222 0.4444444
## [23,] 0.4444444 0.4444444
## [24,] 0.6666667 0.4444444
## [25,] 0.8888889 0.4444444
## [26,] 1.1111111 0.4444444
## [27,] 1.3333333 0.4444444
## [28,] 1.5555556 0.4444444
## [29,] 1.7777778 0.4444444
## [30,] 2.0000000 0.4444444
## [31,] 0.0000000 0.6666667
## [32,] 0.2222222 0.6666667
## [33,] 0.4444444 0.6666667
## [34,] 0.6666667 0.6666667
## [35,] 0.8888889 0.6666667
## [36,] 1.1111111 0.6666667
## [37,] 1.3333333 0.6666667
## [38,] 1.5555556 0.6666667
## [39,] 1.7777778 0.6666667
## [40,] 2.0000000 0.6666667
## [41,] 0.0000000 0.8888889
## [42,] 0.2222222 0.8888889
## [43,] 0.4444444 0.8888889
## [44,] 0.6666667 0.8888889
## [45,] 0.8888889 0.8888889
## [46,] 1.1111111 0.8888889
## [47,] 1.3333333 0.8888889
## [48,] 1.5555556 0.8888889
## [49,] 1.7777778 0.8888889
## [50,] 2.0000000 0.8888889
## [51,] 0.0000000 1.1111111
## [52,] 0.2222222 1.1111111
## [53,] 0.4444444 1.1111111
## [54,] 0.6666667 1.1111111
## [55,] 0.8888889 1.1111111
## [56,] 1.1111111 1.1111111
## [57,] 1.3333333 1.1111111
## [58,] 1.5555556 1.1111111
## [59,] 1.7777778 1.1111111
## [60,] 2.0000000 1.1111111
## [61,] 0.0000000 1.3333333
## [62,] 0.2222222 1.3333333
## [63,] 0.4444444 1.3333333
## [64,] 0.6666667 1.3333333
## [65,] 0.8888889 1.3333333
## [66,] 1.1111111 1.3333333
## [67,] 1.3333333 1.3333333
## [68,] 1.5555556 1.3333333
## [69,] 1.7777778 1.3333333
## [70,] 2.0000000 1.3333333
## [71,] 0.0000000 1.5555556
## [72,] 0.2222222 1.5555556
## [73,] 0.4444444 1.5555556
## [74,] 0.6666667 1.5555556
## [75,] 0.8888889 1.5555556
## [76,] 1.1111111 1.5555556
## [77,] 1.3333333 1.5555556
## [78,] 1.5555556 1.5555556
## [79,] 1.7777778 1.5555556
## [80,] 2.0000000 1.5555556
## [81,] 0.0000000 1.7777778
## [82,] 0.2222222 1.7777778
## [83,] 0.4444444 1.7777778
## [84,] 0.6666667 1.7777778
## [85,] 0.8888889 1.7777778
## [86,] 1.1111111 1.7777778
## [87,] 1.3333333 1.7777778
## [88,] 1.5555556 1.7777778
## [89,] 1.7777778 1.7777778
## [90,] 2.0000000 1.7777778
## [91,] 0.0000000 2.0000000
## [92,] 0.2222222 2.0000000
## [93,] 0.4444444 2.0000000
## [94,] 0.6666667 2.0000000
## [95,] 0.8888889 2.0000000
## [96,] 1.1111111 2.0000000
## [97,] 1.3333333 2.0000000
## [98,] 1.5555556 2.0000000
## [99,] 1.7777778 2.0000000
## [100,] 2.0000000 2.0000000
z <- outer(x, y, cuadri) # intensidad calculada en cada punto de grilla
u<-runif(n*n) # puntos aleatorios x y
graf<-persp(x, y, z,theta = 65,phi=25,zlab="u")
cumplen<-(u<z) # restricción para simular bajo la intensidad
points(trans3d(xx, yy, u,pmat=graf), col = 2+cumplen, pch = 16)
# Simulacion
set.seed(1) # fijo la semilla
cant<-1000 # cant de obs
x<-runif(cant,0,2) # genero x entre 0 y 2
y<-runif(cant,0,2) # genero y entre 0 y 2
u<-runif(cant,0,1) # generu u entre 0 y 1
z<-cuadri(x,y) # evaluo en la funcion de intensidad
cumplen<-(u<z) # Restriccion
sum(cumplen) # cuantos quedan
## [1] 554
# Grafico 100 casos
par(pty="s")
plot(x[cumplen][1:100],y[cumplen][1:100],xlab = "",ylab="",xaxt='n',yaxt='n')
# Grafico todos los casos
par(pty="s")
plot(x[cumplen],y[cumplen],xlab = "",ylab="",xaxt='n',yaxt='n')
Simulacion de Complete Spatial Randomness (CSR) o Independent Random Process (IRP) Funciones de spatstat,core (Uniform) Binomial Point Process
set.seed(1)
cant<-100
sim.BPP<-runifpoint(cant)
class(sim.BPP)
## [1] "ppp"
plot(sim.BPP)
set.seed(2)
sim.BPP<-runifpoint(cant)
plot(sim.BPP)
# muchas realizaciones
sim.BPP<-runifpoint(cant,nsim = 6)
plot(sim.BPP)
class(sim.BPP)
## [1] "ppplist" "solist" "anylist" "listof" "list"
Proceso de poisson homogeneo
set.seed(1)
lambda<-100
sim.CSR<-rpoispp(lambda)
plot(sim.CSR,main = paste("Poisson Homogeneo lambda =",lambda))
set.seed(2)
lambda<-100
sim.CSR<-rpoispp(lambda)
plot(sim.CSR,main = paste("Poisson Homogeneo lambda =",lambda))
Proceso de poisson INhomogeneo
set.seed(1)
sim.INHa <- rpoispp(function(x,y) {800*(x-0.5)^2 + 800*(y-0.5)^2},nsim=1)
plot(sim.INHa)
set.seed(1)
sim.INHb <- rpoispp(function(x,y) {200-800*(x-0.5)^2 + 200-800*(y-0.5)^2},nsim=1)
plot(sim.INHb)
Procesos de Superposición (Poisson)
# a lo bruto
plot(sim.INHa,type="n",main="Superposición")
points(sim.INHa,col="blue")
points(sim.INHb,col="red")
# mas sofisticado
super<-superimpose(sim.INHa,sim.INHb) # superposicion de ppps
plot(super)
Analisis de Concentracion PP Homogeneo
# El proceso CSR
sim.CSR
## Planar point pattern: 91 points
## window: rectangle = [0, 1] x [0, 1] units
plot(sim.CSR)
# Quadrat Count
Q <- quadratcount(sim.CSR, nx= 4, ny=4)
plot(sim.CSR, pch=20, cols="grey70", main=NULL) # Plot points
plot(Q, add=TRUE) # Add quadrat grid
# Quadrat Density
Q.d <- intensity(Q)
plot(intensity(Q, image=TRUE), main=NULL, las=1) # Plot density raster
plot(sim.CSR, pch=1, cex=0.01, col="grey", add=TRUE) # Add points
Test de uniformidad de distribucion PP Homogeneo
QT <- quadrat.test(sim.CSR, nx=4, ny=4)
QT
##
## Chi-squared test of CSR using quadrat counts
##
## data: sim.CSR
## X2 = 16.077, df = 15, p-value = 0.7539
## alternative hypothesis: two.sided
##
## Quadrats: 4 by 4 grid of tiles
plot(QT)
Analisis de Concentracion PP InHomogeneo
# Quadrat Count
Q <- quadratcount(sim.INHa, nx= 4, ny=4)
plot(sim.INHa, pch=20, cols="grey70", main=NULL) # Plot points
plot(Q, add=TRUE) # Add quadrat grid
# Quadrat Density
Q.d <- intensity(Q)
plot(intensity(Q, image=TRUE), main=NULL, las=1) # Plot density raster
plot(sim.INHa, pch=1, cex=0.01, col="grey", add=TRUE) # Add points
Test de uniformidad de distribucion PP InHomogeneo
QT <- quadrat.test(sim.INHa, nx=4, ny=4)
QT
##
## Chi-squared test of CSR using quadrat counts
##
## data: sim.INHa
## X2 = 73.619, df = 15, p-value = 2.009e-09
## alternative hypothesis: two.sided
##
## Quadrats: 4 by 4 grid of tiles
plot(QT)
Proceso de Thomas
par(mfrow=c(2,2))
par(mar=c(1,1,1,1))
set.seed(1)
inten<-30
disper<-0.02
mu<-10
sim.THO <- rThomas(inten,disper,mu)
plot(sim.THO,main=paste("Lambda=",inten,"Mu=",mu,"Sigma=",disper))
set.seed(1)
inten<-30
disper<-0.05
mu<-10
sim.THO <- rThomas(inten,disper,mu)
plot(sim.THO,main=paste("Lambda=",inten,"Mu=",mu,"Sigma=",disper))
set.seed(1)
inten<-10
disper<-0.02
mu<-30
sim.THO <- rThomas(inten,disper,mu)
plot(sim.THO,main=paste("Lambda=",inten,"Mu=",mu,"Sigma=",disper))
set.seed(1)
inten<-10
disper<-0.05
mu<-30
sim.THO <- rThomas(inten,disper,mu)
plot(sim.THO,main=paste("Lambda=",inten,"Mu=",mu,"Sigma=",disper))
par(mfrow=c(1,1))
Proceso de Secuencial Inhibicion Simple (SSI)
par(mfrow=c(2,2))
par(mar=c(1,1,1,1))
set.seed(1)
delta<-0.01
sim.SSI <- rSSI(delta,100)
plot(sim.SSI,main=paste("n=100 delta=",delta))
set.seed(1)
delta<-0.05
sim.SSIa <- rSSI(delta,100)
plot(sim.SSIa,main=paste("n=100 delta=",delta))
set.seed(1)
delta<-0.1
sim.SSI <- rSSI(delta,100)
## Warning in rSSI(delta, 100): Gave up after 1000 attempts with only 69 points
## placed out of 100
plot(sim.SSI,main=paste("n=100 delta=",delta))
set.seed(1)
delta<-0.2
sim.SSIb <- rSSI(delta,100)
## Warning in rSSI(delta, 100): Window is too small to fit 100 points at minimum
## separation 0.2 (absolute maximum number is 31)
## Warning in rSSI(delta, 100): Gave up after 1000 attempts with only 21 points
## placed out of 100
plot(sim.SSIb,main=paste("n=100 delta=",delta))
Procesos de COX
library(RandomFields)
## Loading required package: sp
## Loading required package: RandomFieldsUtils
##
## Attaching package: 'RandomFields'
## The following object is masked from 'package:RandomFieldsUtils':
##
## RFoptions
## The following object is masked from 'package:nlme':
##
## Variogram
library(RandomFieldsUtils)
# Defino el campo medio de la funcion de intensidad
m <- as.im(function(x, y){5 - 1.5 * (x - 0.5)^2 + 2 * (y - 0.5)^2}, W=owin())
class(m)
## [1] "im"
plot(m)
# Genero un Proceso de Cox Log Gausiano a partir del campo medio generado previamente
X <- rLGCP("gauss", m, var=0.15, scale =0.5)
# Grafico del campo aleatorio + la realizacion del proceso
plot(attr(X, "Lambda"))
points(X)
# Probemos tra ventana
#
ho <- owin(poly=list(list(x=c(0,2,1,0), y=c(0,0.3,1,0.8)),
list(x=c(0.6,0.4,0.4,0.6), y=c(0.2,0.2,0.4,0.4))))
m <- as.im(function(x, y){5 - 1.5 * (x - 0.5)^2 + 2 * (y - 0.5)^2}, W=ho)
plot(m)
La Funcion G y F en CSR
par(mfrow=c(1,3))
plot(sim.CSR)
env.CSR <- envelope(sim.CSR, fun = Gest, nsim = 200,nrank=5)
## Generating 200 simulations of CSR ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198. 200.
##
## Done.
plot(env.CSR)
env.CSR <- envelope(sim.CSR, fun = Fest, nsim = 200,nrank=5)
## Generating 200 simulations of CSR ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198. 200.
##
## Done.
plot(env.CSR)
par(mfrow=c(1,1))
La Funcion G y F en Thomas (Concentracion)
par(mfrow=c(1,3))
plot(sim.THO)
env.THO <- envelope(sim.THO, fun = Gest, nsim = 200,nrank=5)
## Generating 200 simulations of CSR ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198. 200.
##
## Done.
plot(env.THO)
env.THO <- envelope(sim.THO, fun = Fest, nsim = 200,nrank=5)
## Generating 200 simulations of CSR ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198. 200.
##
## Done.
plot(env.THO)
par(mfrow=c(1,1))
La Funcion G y F en SSI (Regularidad)
par(mfrow=c(1,3))
sim.SSI<-sim.SSIa
plot(sim.SSI)
env.SSI <- envelope(sim.SSI, fun = Gest, nsim = 200,nrank=5)
## Generating 200 simulations of CSR ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198. 200.
##
## Done.
plot(env.SSI)
env.SSI <- envelope(sim.SSI, fun = Fest, nsim = 200,nrank=5)
## Generating 200 simulations of CSR ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198. 200.
##
## Done.
plot(env.SSI)
par(mfrow=c(1,1))
La Funcion G y F en CSR
par(mfrow=c(1,3))
plot(sim.CSR)
env.CSR <- envelope(sim.CSR, fun = Kest, nsim = 200,nrank=5)
## Generating 200 simulations of CSR ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198. 200.
##
## Done.
plot(env.CSR)
env.CSR <- envelope(sim.CSR, fun = Lest, nsim = 200,nrank=5)
## Generating 200 simulations of CSR ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198. 200.
##
## Done.
plot(env.CSR, . -r ~ r)
par(mfrow=c(1,1))
La Funcion K y L en Thomas
par(mfrow=c(1,3))
plot(sim.THO)
env.THO <- envelope(sim.THO, fun = Kest, nsim = 200,nrank=5)
## Generating 200 simulations of CSR ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198. 200.
##
## Done.
plot(env.THO)
env.THO <- envelope(sim.THO, fun = Lest, nsim = 200,nrank=5)
## Generating 200 simulations of CSR ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198. 200.
##
## Done.
plot(env.THO, . -r ~ r)
par(mfrow=c(1,1))
La Funcion K y L en SSI
par(mfrow=c(1,3))
sim.SSI<-sim.SSIa
plot(sim.SSI)
env.SSI <- envelope(sim.SSI, fun = Kest, nsim = 200,nrank=5)
## Generating 200 simulations of CSR ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198. 200.
##
## Done.
plot(env.SSI)
env.SSI <- envelope(sim.SSI, fun = Lest, nsim = 200,nrank=5)
## Generating 200 simulations of CSR ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198. 200.
##
## Done.
plot(env.SSI, . -r ~ r)
par(mfrow=c(1,1))
Fry Plots o graficos de Petterson
par(mfrow=c(2,3))
par(mar=c(1,1,1,1))
plot(sim.CSR)
plot(sim.THO)
plot(sim.SSIb)
fryplot(sim.CSR)
fryplot(sim.THO)
fryplot(sim.SSIb)
par(mfrow=c(1,1))
Estimcion de la intensidad para poisson con modelo parametrico (funcion ppm)
Caso Homogeneo
# modelo homogeneo
ajus.CSR.1<-ppm(sim.CSR~1)
ajus.CSR.1
## Stationary Poisson process
## Intensity: 91
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## log(lambda) 4.51086 0.1048285 4.305399 4.71632 *** 43.03086
# modelo cinhomogeneo on tendencia espacial lineal
ajus.CSR.2<-ppm(sim.CSR~x+y)
ajus.CSR.2
## Nonstationary Poisson process
##
## Log intensity: ~x + y
##
## Fitted trend coefficients:
## (Intercept) x y
## 4.52381507 0.06116668 -0.08803713
##
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## (Intercept) 4.52381507 0.2769819 3.9809404 5.0666897 *** 16.3325273
## x 0.06116668 0.3633154 -0.6509185 0.7732518 0.1683569
## y -0.08803713 0.3633713 -0.8002318 0.6241576 -0.2422787
plot(ajus.CSR.2,pause=FALSE,superimpose = FALSE)
plot(ajus.CSR.2,pause=FALSE,how="persp", theta=-30,phi=40,d=4)
plot(ajus.CSR.2,pause=FALSE,how="contour", theta=-30,phi=40,d=4)
pred2<-predict(ajus.CSR.2)
plot(pred2)
Caso InHomogeneo
# modelo homogeneo
ajus.INH.1<-ppm(sim.INHa~1)
ajus.INH.1
## Stationary Poisson process
## Intensity: 126
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## log(lambda) 4.836282 0.08908708 4.661674 5.010889 *** 54.28713
# modelo cinhomogeneo on tendencia espacial lineal
ajus.INH.2<-ppm(sim.INHa~x+y)
ajus.INH.2
## Nonstationary Poisson process
##
## Log intensity: ~x + y
##
## Fitted trend coefficients:
## (Intercept) x y
## 4.68074523 -0.08171004 0.38018558
##
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## (Intercept) 4.68074523 0.2413345 4.2077383 5.1537521 *** 19.3952605
## x -0.08171004 0.3087647 -0.6868778 0.5234577 -0.2646353
## y 0.38018558 0.3097641 -0.2269410 0.9873121 1.2273389
plot(ajus.INH.2,pause=FALSE,superimpose = FALSE)
plot(ajus.INH.2,pause=FALSE,how="persp", theta=-30,phi=40,d=4)
plot(ajus.INH.2,pause=FALSE,how="contour", theta=-30,phi=40,d=4)
pred2<-predict(ajus.INH.2)
plot(pred2)
# modelo cinhomogeneo on tendencia espacial NO lineal
require(splines)
## Loading required package: splines
#ajus.INH.3<-ppm(sim.INHa~bs(x,df=2)+bs(y,df=2),Poisson())
ajus.INH.3<-ppm(sim.INHa~poly(x,2)+poly(y,2),Poisson())
ajus.INH.3
## Nonstationary Poisson process
##
## Log intensity: ~poly(x, 2) + poly(y, 2)
##
## Fitted trend coefficients:
## (Intercept) poly(x, 2)1 poly(x, 2)2 poly(y, 2)1 poly(y, 2)2
## 4.6815325 -0.3972158 15.8736932 2.8181625 14.0666392
##
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## (Intercept) 4.6815325 0.1019326 4.481748 4.881317 *** 45.9277186
## poly(x, 2)1 -0.3972158 2.5786078 -5.451194 4.656763 -0.1540427
## poly(x, 2)2 15.8736932 2.8529773 10.281960 21.465426 *** 5.5639046
## poly(y, 2)1 2.8181625 2.6344315 -2.345228 7.981553 1.0697422
## poly(y, 2)2 14.0666392 2.8831711 8.415728 19.717551 *** 4.8788777
plot(ajus.INH.3,pause=FALSE,superimpose = FALSE)
plot(ajus.INH.3,pause=FALSE,how="persp", theta=-30,phi=40,d=4)
plot(ajus.INH.3,pause=FALSE,how="contour", theta=-30,phi=40,d=4)
pred3<-predict(ajus.INH.3)
plot(pred3)
Como funciona KDE
library(manipulate)
x<-c(10,13,14)
GrafKDE<-function(h)
{
plot(density(x,bw=h,from=8,to=16),main=paste("Ventana =",h),ylab="Intensidad",ylim=c(0,0.8),lwd=3,xlab="")
points(x,c(0,0,0),pch=16,col="green")
n1<-function(x){dnorm(x,10,h)/3}
n2<-function(x){dnorm(x,13,h)/3}
n3<-function(x){dnorm(x,14,h)/3}
curve(n1, 7, 17, lwd=2, axes = FALSE, xlab = "", ylab = "",add=T,col="green",lty=3)
curve(n2, 7, 17, lwd=2, axes = FALSE, xlab = "", ylab = "",add=T,col="green",lty=3)
curve(n3, 7, 17, lwd=2, axes = FALSE, xlab = "", ylab = "",add=T,col="green",lty=3)
}
GrafKDE(0.5)
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "axes" is not a graphical
## parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "axes" is not a graphical
## parameter
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "axes" is not a graphical
## parameter
#manipulate(GrafKDE(h=h),h=slider(0.2,2,step=0.1))
Aplicacion de KDE a Poisson Homogeneo
require(splancs)
## Loading required package: splancs
##
## Spatial Point Pattern Analysis Code in S-Plus
##
## Version 2 - Spatial and Space-Time analysis
# eleccion de ventana por diggle
mserw <- bw.diggle(sim.CSR)
bw <- as.numeric(mserw)*3
#
den.CSR<-density(sim.CSR,bw)
plot(den.CSR)
points(sim.CSR)
Aplicacion de KDE a Poisson InHomogeneo
require(splancs)
# eleccion de ventana por diggle
mserw <- bw.diggle(sim.INHa)
bw <- as.numeric(mserw)*3
#
den.INH<-density(sim.INHa,bw)
plot(den.INH)
points(sim.INHa)
Elegiendo la ventana optima para Poisson Inhomogeneo
require(spatstat.utils)
## Loading required package: spatstat.utils
#datos<-sim.CSR
datos<-sim.INHa
poligono<-matrix(c(0,0,0,1,1,0,1,1,0,0),4,2,byrow = T)
## Warning in matrix(c(0, 0, 0, 1, 1, 0, 1, 1, 0, 0), 4, 2, byrow = T): data length
## [10] is not a sub-multiple or multiple of the number of rows [4]
Mse2d <- mse2d(as.points(datos), poligono, nsmse=50, range=1)
#png("/home/andresfaral/Dropbox/Estadistica Espacial/KDE INH CV")
plot(Mse2d$h,Mse2d$mse, type="l",xlab="Ventana (h)", ylab="Error Estimado")
arrows(0.15,9,0.15,-1.5,lty=1,col="green",lwd=3)
#dev.off()
bw<-min(Mse2d$h)
bw<-0.03
den.INH<-density(datos,bw)
png("/home/andresfaral/Dropbox/Estadistica Espacial/KDE INH chico")
plot(den.INH)
points(datos)
dev.off()
## png
## 2
Elegiendo la ventana optima para Thomas
datos<-sim.THO
poligono<-matrix(c(0,0,0,1,1,0,1,1,0,0),4,2,byrow = T)
## Warning in matrix(c(0, 0, 0, 1, 1, 0, 1, 1, 0, 0), 4, 2, byrow = T): data length
## [10] is not a sub-multiple or multiple of the number of rows [4]
Mse2d <- mse2d(as.points(datos), poligono, nsmse=50, range=1)
bw<-Mse2d$h[which.min(Mse2d$mse)]
png("/home/andresfaral/Dropbox/Estadistica Espacial/KDE THO CV")
plot(Mse2d$h,Mse2d$mse, type="l",xlab="Ventana (h)", ylab="Error Estimado")
arrows(bw,-1,bw,-4.5,lty=1,col="green",lwd=2)
dev.off()
## png
## 2
den.THO<-density(datos,bw)
png("/home/andresfaral/Dropbox/Estadistica Espacial/KDE THO")
plot(den.THO)
points(datos)
dev.off()
## png
## 2
Procesos Marcados
library(gstat)
##
## Attaching package: 'gstat'
## The following object is masked from 'package:spatstat.core':
##
## idw
library(lattice)
##
## Attaching package: 'lattice'
## The following object is masked from 'package:spatstat.core':
##
## panel.histogram
library(sp)
data(meuse)
class(meuse)
## [1] "data.frame"
coordinates(meuse) <- c("x", "y")
# EDA
spplot(meuse, "zinc", do.log = T, colorkey = TRUE)
bubble(meuse, "zinc", do.log = T, key.space = "bottom")
# fitting
xyplot(log(zinc) ~ sqrt(dist), as.data.frame(meuse))
zn.lm <- lm(log(zinc) ~ sqrt(dist), meuse)
meuse$fitted.s <- predict(zn.lm, meuse) - mean(predict(zn.lm,
meuse))
meuse$residuals <- residuals(zn.lm)
spplot(meuse, c("fitted.s", "residuals"))
Interpolacion por IDW
data(meuse.grid)
class(meuse.grid)
## [1] "data.frame"
coordinates(meuse.grid) <- c("x", "y")
meuse.grid <- as(meuse.grid, "SpatialPixelsDataFrame")
idw.out <- gstat::idw(zinc ~ 1, meuse, meuse.grid, idp = 0.5)
## [inverse distance weighted interpolation]
plot(idw.out)
Interpolacion con Modelado Lineal
#zn.lm <- lm(log(zinc) ~ sqrt(dist), meuse)
#zn.lm <- lm(zinc ~ sqrt(dist), meuse)
zn.lm <- lm(zinc ~ poly(x, y, degree = 2), meuse)
meuse.grid$pred <- predict(zn.lm, meuse.grid)
meuse.grid$se.fit <- predict(zn.lm, meuse.grid, se.fit = TRUE)$se.fit
plot(meuse.grid["pred"])
#image(meuse.grid["pred"])
El Variograma de zinc
hist(dist(coordinates(meuse)))
hscat(zinc ~ 1, meuse, seq(0,2000,length.out = 10))
#
sel <- plot(variogram(zinc ~ 1, meuse, cloud = FALSE))
sel
El Variograma de residuos
hscat(residuals ~ 1, meuse, seq(0,2000,length.out = 10))
#
sel <- plot(variogram(residuals ~ 1, meuse, cloud = FALSE))
sel
#
sel2<-variogram(log(zinc) ~ sqrt(dist), meuse)
plot(sel2)
#
zn.lm2 <- lm(zinc ~ dist+poly(x, y, degree = 2), meuse)
meuse$fitted2 <- predict(zn.lm2, meuse) - mean(predict(zn.lm2,
meuse))
meuse$residuals2 <- residuals(zn.lm2)
hscat(residuals2 ~ 1, meuse, seq(0,2000,length.out = 10))
#
sel <- plot(variogram(residuals2 ~ 1, meuse, cloud = FALSE))
sel
Procesos Marcados
SIN Autocorrelacion
require(gstat)
#
set.seed(1)
sim.CSR.M<-sim.CSR
X<-rnorm(sim.CSR.M$n)
X<-as.numeric(scale(X))
marks(sim.CSR.M)<-X
sim.CSR.M
## Marked planar point pattern: 91 points
## marks are numeric, of storage type 'double'
## window: rectangle = [0, 1] x [0, 1] units
sim.CSR.M.df<-as.data.frame(sim.CSR.M)
sim.CSR.M.spdf<-sim.CSR.M.df
coordinates(sim.CSR.M.spdf) <- c("x", "y") # como spatialpoints dataframe
plot(sim.CSR.M)
spplot(sim.CSR.M.spdf, colorkey = TRUE)
### Vaiograma
hscat(marks ~ 1, sim.CSR.M.spdf, seq(0,1,length.out = 11))
variog<-variogram(marks ~ 1, sim.CSR.M.spdf)
plot(variog,type="l")
CON Autocorrelacion
A proṕosito NO uso ni gstat ni fields
Genero valores al azar y luego hago muchas pasadas (pp) en las que elijo un punto al azar y modifico su marca sumandole una fracción (ff) de los valores de los (Kk) puntos cercanos
require(nabor) # Libreria para calcular vecinos mas cercanos
## Loading required package: nabor
sim.CSR.M2<-sim.CSR.M
n<-sim.CSR.M2$n
pp<-500
ff<-0.5
kk<-3
#
coordenadas<-coords(sim.CSR.M2)
vecinos<-knn(coordenadas,coordenadas,k=kk+1)$nn.idx[,2:(kk+1)] # busca vecinos
X.sac<-marks(sim.CSR.M2)
for (i in 1:pp)
{
cual<-sample(1:n,1)
X.sac[cual]<-X.sac[cual]+ff*mean(X.sac[vecinos[cual,]])
}
X.sac<-as.numeric(scale(X.sac))
marks(sim.CSR.M2)<-X.sac
sim.CSR.M2.df<-as.data.frame(sim.CSR.M2)
sim.CSR.M2.spdf<-sim.CSR.M2.df
coordinates(sim.CSR.M2.spdf) <- c("x", "y") # como spatialpoints dataframe
plot(sim.CSR.M2)
spplot(sim.CSR.M2.spdf, colorkey = TRUE)
### Vaiograma
hscat(marks ~ 1, sim.CSR.M2.spdf, seq(0,1,length.out = 11))
variog<-variogram(marks ~ 1, sim.CSR.M2.spdf)
plot(variog,type="l")
Indice de Moran
require(ape)
## Loading required package: ape
##
## Attaching package: 'ape'
## The following object is masked from 'package:splancs':
##
## zoom
## The following objects are masked from 'package:spatstat.geom':
##
## edges, rotate
#
datos<-sim.CSR.M
datos.dists <- as.matrix(dist(cbind(sim.CSR.M$x, sim.CSR.M$y)))
datos.dists.inv <- 1/datos.dists
diag(datos.dists.inv) <- 0
Moran.I(marks(datos),datos.dists.inv)
## $observed
## [1] -0.04942335
##
## $expected
## [1] -0.01111111
##
## $sd
## [1] 0.01848804
##
## $p.value
## [1] 0.03824008
#
datos<-sim.CSR.M2
datos.dists <- as.matrix(dist(cbind(sim.CSR.M$x, sim.CSR.M$y)))
datos.dists.inv <- 1/datos.dists
diag(datos.dists.inv) <- 0
Moran.I(marks(datos),datos.dists.inv)
## $observed
## [1] 0.1257112
##
## $expected
## [1] -0.01111111
##
## $sd
## [1] 0.01842769
##
## $p.value
## [1] 1.130207e-13
Kriging
require(gstst)
## Loading required package: gstst
## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'gstst'
data(meuse.grid)
class(meuse)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
class(meuse.grid)
## [1] "data.frame"
coordinates(meuse.grid) <- c("x", "y")
kri.meuse<- krige(zinc~1, meuse, meuse.grid)
## [inverse distance weighted interpolation]
class(meuse.grid)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
plot(kri.meuse)
spplot(kri.meuse["var1.pred"])
#
M2.grid<-expand.grid(x=seq(0,1,length.out = 100),y=seq(0,1,length.out = 100))
coordinates(M2.grid) <- c("x", "y")
class(M2.grid)
## [1] "SpatialPoints"
## attr(,"package")
## [1] "sp"
m <- vgm(.59, "Sph", 874, .04)
kri.M2 <- krige(marks ~ 1, sim.CSR.M2.spdf,M2.grid,model=m)
## [using ordinary kriging]
plot(kri.M2)
plot(sim.CSR.M2)
spplot(kri.M2["var1.pred"], colorkey = TRUE)